function [] = gendata(hyp_mean)
% a script to generate data by a Gaussian process directly
% by Mark Norrish, 2011
% hyp_mean = [0;0];
range = 5; % i.e. a grid from 1 to range
n = range*range;
n_class = 3;
nfactor = 4; % the higher this parameter, the easier the problem

D1 = repmat((1:range)', range, 1)/nfactor;
D2 = kron(1:range,ones(1, range))'/nfactor;
D = [D1 D2];
%D = (1:range)' / nfactor;
func = @covSEiso; % use covSEiso

if nargin == 0
  hyp_mean = zeros(str2num(func()),1);
end

%hyp_mean = [0 0 0]';
%hypmean=[0, 1:x,0]
hyp_var = [1 1]';
%hyps = normrnd(repmat(hyp_mean, 1, n_class), repmat(hyp_var, 1, n_class)); % randomly generate hyps
hyps = zeros(str2num(func()), n_class);
hyps(:,1) = hyp_mean;

sigma = 1e-4;
K = zeros(n,n,n_class);

for i = 1:n_class
  K(:,:,i) = func(hyps(:,i), D, D) + sigma*eye(n);
end

cholk = zeros(n*n_class);
for c = 1:n_class
  cholk(1+(c-1)*n:c*n,1+(c-1)*n:c*n) = chol(K(:,:,c));
end

F = cholk'*normrnd(zeros(n*n_class,1), ones(n*n_class,1)); % because mu is zero

Fdash = reshape(F, n, n_class)';

Y = (Fdash==repmat(max(Fdash),n_class,1));
Y = sum(Y.*repmat((1:n_class)',1,n));

if ~all(ismember(1:n_class, Y))
  gendata(hyp_mean);
else
  dlmwrite('gp_data', [D Y']);
  dlmwrite('gp_function', F);
  dlmwrite('gp_hyps', hyps);
end